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Abstract. The collective behavior of a two-dimensional wet granular cluster under 
horizontal swirling motions is investigated experimentally. Depending on the balance 
between the energy injection and dissipation, the cluster evolves into various 
nonequilibrium stationary states with strong internal structure fluctuations with time. 
Quantitative characterizations of the fluctuations with the bond orientational order 
parameter qe reveal power spectra of the form /“ with the exponent a closely related 
to the stationary states of the system. In particular, 1 // type of noise with a«-l 
emerges as melting starts from the free surface of the cluster, suggesting the possibility 
of using 1// noise as an indicator for phase transitions in systems driven far from 
thermodynamic equilibrium. 


1. Introduction 

1 // noise (also called flicker or pink noise) is a class of signal exhibiting power spectrum 
/“ with exponent a ~ —1, typically in the range —1.4 < a < —0.8 p]. Due to its 
ubiquity in nature, 1// noise has drawn much more attentions from various communities 
for almost a century, in comparison to white (ct = 0), Brownian (ct = —2) as well as 
diffusion noises [2j |3j IU E] . The first and most extensively investigated noise of this class 
is electric noise in conducting or semiconducting materials P El El IB]. Later on, the 
ubiquity of 1 // noise is unveiled as a surprisingly large number of systems are found to 
fall into this class, ranging from blinking of stars and sunspot activity [2] in astrophysics, 
earthquake triggering [9] and undersea ocean currents in geophysics [TO], to the loudness 
fluctuations of music EH, gene expression [12], and human cognition pam. 

Because of this ubiquity, it is important to understand why such a common feature 
persists and to explore whether there is a universal mechanism behind or not. Following 
the widely used the Bernamont-Surdin-McWhorter (BSM) model jl5l fTHl fl7l fl8] . 1/f 
noise can be considered as a superposition of independent events with Lorentzian spectra 
and a broad distribution of relaxation time. Although it was originally introduced to 
explain electric noises, this model has been successfully applied to a large variety of 
systems exhibiting 1/f noise [U[19]. Nevertheless, recent investigations also reveal that 
the assumption of elementary events with Lorentzian spectra may not be necessary, since 
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models based on random or point process [221 21] have also been reported to produce 
1 // noise. 

In order to provide a general understanding of 1// noise, Bak et al. proposed the 
concept of self-organized criticality (SOC) j22j 23], claiming that in spatially extended 
dissipative systems 1/f noise can be viewed as an indication of self-organized critical 
state. In contrast to the critical state in equilibrium thermodynamics, the system is 
self-organized, i.e. no external tuning is needed. Using a cellular automata (CA) model 
describing the avalanches in a pile of sand, they demonstrated that the scale invariance 
in time is associated with fractal structures, i.e. scale invariance in space. Triggered by 
this concept, a tremendous amount of investigations have been conducted to test the 
universality of SOC [T9], 2T 5 [25] . However, the CA model proposed by BTW was soon 
found to produce 1/f 2 noise in both one- and two-dimensions [26], 27] and thus cannot 
directly be used to explain 1/f noise. Moreover, contradictory results were found in 
model experiments on the avalanches of a sandpile [231128], [29]. Those investigations led 
to the conclusion that SOC is not as universal as it was claimed. Instead, its validity 
relies on the dissipation mechanisms of specific systems [30], and also on the way of 
analysis [an US]. Thus from the SOC perspective, the ubiquitous 1/f noise is still not 
completely understood and whether a general theory for 1/f noise exists or not is still 
unclear. Nevertheless, the advance on the relation between the invariance of time and 
space gives us the hint that noise may provide useful insights into the self-organized 
stationary states of systems driven far from thermodynamic equilibrium [T9] . 

Wet granular matter, such as wet sand used to build sand sculptures, has been 
drawing more and more interest from both physical and engineering communities in 
the past decades [33] EH [35j. This is not only because of its ubiquity in nature, 
industries and our daily lives, but also due to the fact that it can be treated as a 
nonequilibrium model system with granular particles replacing molecules and liquid 
bridges formed between adjacent particles replacing molecular bonds. Because of the 
strong energy dissipation associated with wet impacts [36], continuous energy injection is 
necessary to maintain a certain stationary state of such a nonequilibrium system. Former 
investigations have revealed that the dissipative energy scale arising from the rupture 
of liquid bridges between adjacent particles plays an important role in determining the 
collective behavior of wet granular matter, such as melting [37] EE], clustering [39], 
and phase separation [TO] SI]- More recently, an analogue of surface melting was 
found in agitated wet granular matter mm, suggesting that existing knowledge on 
phase transitions in equilibrium systems can be extended to explain the wide-spreading 
nonequilibrium systems in nature. Here, focusing on the noise spectra extracted the 
nonequilibrium stationary states of a wet granular cluster, I demonstrate the possibility 
of using 1/f type noise as an indicator for phase transitions far from thermodynamic 
equilibrium. 
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Figure 1 . A sketch of the experimental setup. The swirling table provides circular 
motion of the whole container in a plane perpendicular to the gravitational acceleration 
g. Both the CCD camera and the LED lamp are mounted in the co-moving frame of 
the swirling table, (a) A typical snapshot of a wet granular “crystal” with defects. 
The circular arrow denotes the swirling motion of the container, (b) The processed 
image with particles color coded according to their local structures: Red, green and 
blue represent hexagonal, square and chain structures, respectively. 


2. Methods 

The experimental setup is sketched in Fig. CD Various number N of polished soda 
lime glass beads (SiLibeads P) with a density of p p = 2.58 g/cm 3 and a diameter of 
d — 4 ± 0.02 mm are used as granular sample. After being mixed with a certain volume 
V\ of purified water (specific resistance 18.2 MDcm, LaborStar TWF), the particles 
are added into a cylindrical polytetrafluoroethylcne (PTFE) container with an inner 
diameter of D = 102 mm and a height of 6 mm. The latter confinement ensures a 
strict mono-layer of particles. The liquid content W = V\/V s with V s the total volume of 
the spheres is fixed at 0.02, which ensures short ranged attractive interactions between 
adjacent particles via the formation of liquid bridges [Hj. The container is closed 
with a glass lid to keep the liquid content constant. During the experiments, the 
glass lid is slightly heated to avoid water vapor condensed there. The wet granular 
sample is illuminated from below with a LED array, and viewed from top with a 
CCD camera (Lumenera LU125M) mounted in the comoving frame [see Fig. Q3(a) for 
a typical snapshot]. The swirling table is leveled within 5.7 x 10” 3 degrees to avoid 
the influence from gravity. The horizontal swirling motion, i.e. a superposition of two 
perpendicular sinusoidal vibrations with identical amplitude A and a phase difference 
of 7t/2, is provided by a mechanical orbit shaker (Thermolyne, AROS160). The amount 
of energy injection into the system is adjusted by varying the driving frequency /a at a 
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fixed amplitude A = 31.8 mm. With a computer controlled precision resistance decay 
box (Burster 1424), the swirling frequency is controlled with an accuracy of ~ 10~ 4 Hz. 
fd is measured by tracing a fixed particle on the swirling table with another CCD camera 
(Lumenera LU075M) mounted in the lab frame. A computer program is developed to 
adjust the swirling frequency in steps and capture images with a fixed rate of 5 frames 
per second. 

Using an image processing procedure based on the Hough transformation [45], 1 
determine the positions of all particles in each frame captured. For each particle, the 
bond orientational order parameters (BOOP) are calculated with [46, 471 HE] 
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where Q\ m = (Qim(A)) is an average of the local order parameter <5i m (r) = Y\ m (6(r), f>(r)) 
over all bonds connecting this particle to its nearest neighbors, which are identified with 
a critical bond length r b = 1.35 d. Here Yi m (9(r), <f{f)) corresponds to the spherical 
harmonics of a bond located at r. Finally, 1 determine the local structure of each 
particle by comparing the order parameter q e with the standard values for perfectly 
hexagonal, square, line structures, as well as for isolated particles. Here, q$ is chosen 
as the order parameter because of its sensitivity to the hexagonal order. As shown in 
Fig. m, 4 ie particles are colored by their local structures after the above analysis in 
order to visualize the local structure changes. 

The collective behavior of a wet granular cluster is obtained through initializing the 
system at a relatively large fd = 1.745 Hz for 2 hours followed by decreasing fd with a 
step of A fd = 0.062 Hz and fixed waiting time At = 700 s. The long initial time ensures 
a homogeneous wetting liquid distribution and a stable temperature gradient to avoid 
the condensation of water on the lid. Immediately after the lowest driving frequency is 
reached, fd is swept up with the same A/ d and At. The whole cycle is repeated seven 
times for each particle number N. 


3. Nonequilibrium stationary states 

Figure [2] shows the collective behavior of a wet granular cluster at four different driving 
frequencies and corresponding internal structure fluctuations. The internal structure is 
characterized by Pq, the percentage of particles in a local hexagonal structure. 

At fd = 1.161 Hz, the cluster is composed of a few “crystals” loosely connected with 
each other. Here, “crystals” refer to clusters composed of particles in a local hexagonal 
packing. As the positions of the cluster center [see the cloud of dots in Fig. [2](a)] indicate, 
the mobility of the cluster is low. This suggests that the external driving force is not 
sufficient for the “crystals” to overcome the friction from the ground and move freely. 
The P 6 order parameter stays at a high value close to 0.9, indicating that almost all 
particles are kept in a hexagonal structure. Thus, this state is called initial crystalline 
state. The small fluctuations of P 6 mainly arise from the interactions between the loosely 
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Figure 2. Sample images of a wet granular cluster in various nonequilibrium 
stationary states and corresponding internal structure fluctuations with time. Upper 
panels: Snapshots of the cluster in its initial crystalline (a), amorphous (b), surface 
melting (c), and completely melted (d) states with the positions of its center of mass 
(c.m.) in the co-moving frame over the whole observation period overlaid as a cloud 
of orange dots. Middle panels: Corresponding internal structures of the cluster based 
on an analysis with BOOP. The color code is the same as described in the caption 
of Fig. [1] Isolated particles are colored in black. Lower panel: The fluctuations of 
the order parameter Pg, i.e. fraction of particles in the hexagonal structure (in red), 
with time. The total number of particles N = 406, corresponding to an area fraction 
$ = 0.624. 

connected “crystals” and the frequent attaching and detaching of individual particles 
from them. 

As the driving frequency / c j increases to 1.285 Hz, the enhanced mobility of the 
particles leads to the merging of the “crystals” into a single one with dislocation defects. 
The defect lines, along which local structures of particles deviate from hexagons to 
squares, suggest a certain amount of potential energy being stored there. This energy 
is obtained from the frequent collisions between the “crystal” and the container wall, 
because of the higher mobility of the whole cluster in comparison to the initial crystalline 
state. Note that the state with defects is unstable: A small perturbation will lead to a 
healing of the defects due to the cohesion between adjacent particles. As the ground state 
with a perfect hexagonal packing is not easily achieved, the healing process typically 
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ends up with a metastable state corresponding to a local potential energy minimum. 
The corresponding fluctuations of P 6 are much stronger than the initially crystalline 
state, because the defects associated structure changes occur at the length scale of the 
whole “crystal”. Moreover, the fluctuations occur at time scales much longer than the 
swirling period, i.e. l//d- 

Such type of fluctuations persist as the driving frequency increases, as long as the 
cluster stays in the amorphous state. As /d increases further to 1.409 Hz, an abrupt 
transition into the surface melting state arises, owing to a balance between the energy 
injection and the rupture energy of a liquid bridge [43j. This transition is manifested 
by the collapse of the voids inside the “crystal” together with the emerging of a liquid 
like layer covering the crystalline core, as shown in Fig. [2](c) and (g). Consequently, the 
whole cluster moves in a circular trajectory with permanent contact to the container 
wall (Note that the shortest distance from the circle to to the container wall matches 
the mean radius of the cluster). This feature indicates that the hard impacts between 
the cluster and the container wall with a finite effective coefficient of restitution (COR) 
is now replaced with soft ones (i.e. zero COR), because the liquid like layer acts as an 
effective damper. 

It is remarkable to see that the fluctuations of P 6 also behave dramatically different 
in the surface melting state in comparison to the amorphous one: The large fluctuations 
at long time scales in the amorphous state are not present any more. This can be 
understood from the fact that the energy injection and dissipation in this state is 
localized in the molten surface layer of the cluster, because of the much more frequent 
collisions there. Hence, the metastability associated with the change of the crystalline 
structures at a length scale of the whole cluster is suppressed, leading to a white noise 
type fluctuations. Further increasing /d leads to an growth of the molten layer thickness, 
until the whole cluster behaves like a liquid at /a = 1.719 Hz. Along with the growth of 
the molten layer thickness, the center of the more deformable cluster is driven closer to 
the rim of the container, as a comparison between Fig. [2](c) and (d) reveals. Moreover, 
the transition into the liquidlike state leads to a decrease of the average Pq, although 
the behavior of the fluctuations persists. 

The various stationary states and the associated structure fluctuations are robust: 
Variations of A/d, At by at least one order of magnitude and a change the sweeping 
direction of /d yield the same behavior. The lower limit of /d is selected such that the 
wet granular assemblies are immobile, i.e. the driving force cannot overcome the total 
frictional force from the ground. 

4. Power spectrum of internal structure fluctuations 

From the perspective of 1// noise, it is intuitive to analyze the power spectrum of the 
internal structure fluctuations in various nonequilibrium stationary states of the wet 
granular cluster. As shown in Fig. [3l the power spectra of P 6 exhibit predominately 
power law behavior for a wide range of /d. In the initial crystalline or amorphous state 
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Figure 3. Normalized power spectrum of Pq at various driving frequency fd- The 
measurement is performed by decreasing /<j. Each power spectrum is an average over 7 
experimental runs. To enhance visibility, the data are smoothed with running average. 
The dash lines (corresponding to power law exponent —0.3 and —1.5) are drawn to 
guide the eyes. Other parameters are the same as described in the caption of Fig. [2] 


(/d < 1.285 Hz), all the power spectra decay with an exponent of a ~ —1.5. The peak 
at / = /a indicates the influence from the driving frequency, as the energy injection is 
associated with periodic collisions between the cluster and the container wall. The peak 
becomes less and less pronounced as /d grows, which can be attributed to the “softening” 
of the impact: As the “crystal” becomes more amorphous, the effective COR for its 
impact with the container wall decreases, and the internal structure change of the cluster 
becomes more susceptible. This trend will eventually lead to more random collisions with 
the container wall and consequently a vanishing influence from the container wall. As 
a = —1.5 is typically associated with noise generated by diffusion mechanisms MS, 
it is interesting to explore possible diffusion mechanisms of driven wet granular particles 
in an amorphous state. For a monolayer of noncohesive dry granular particles, former 
investigations (see e.g. Reis et al. BSD reveal a non-trivial subdiffusive behavior in 
the vicinity of melting transition, which arises from the caging dynamics. However, 
it is still unclear what the diffusion mechanisms of wet granular particles in various 
nonequilibrium stationary states are and how it is related to the power spectra observed 
here. Further clarifications of these questions require a detailed characterization of the 
particle mobility and will be a focus of further investigations. 

The exponent a decreases to ~ —1.2 as /d increases to 1.347 Hz, suggesting a 1 // 
type noise in the vicinity of melting transition. In the surface melting regime with even 
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Figure 4. (a) Power law exponent a as a function of the driving frequency /d for 
area fraction $ = 0.624 (N = 406). Upper and lower triangles represent data obtained 
by increasing and decreasing /d, respectively, (b) a as a function of /d for various 4>. 
Solid lines are fits to the data for determining the threshold frequency / t h (see text for 
details). Inset shows the threshold frequency /th as a function of $. 


higher /a, a grows continuously until it saturates at ~ —0.3, suggesting a trend toward 
a white noise type of fluctuations without any long time correlations. Meanwhile, a 
cut-off frequency of 0.1 Hz, below which the power spectrum decays much slower, 
arises. 

As the power law exponent a is closely related to the stationary states of the system, 
it is intuitive to use it as an order parameter for the melting transition. Figure|U(a) (lower 
triangle) shows the fitted power law exponent as a function of the driving frequency for 
the spectrum shown in Fig. [3j The transition from an amorphous state to a melted 
state can be clearly distinguished from the deviation of a from its initial low value. 
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To avoid the influence from the peak at / d , I fit the spectrum with an upper limit 
/lim, which is determined by decreasing / lim in steps from / d until the standard error 
of the fit converges. Because of the logarithmic scale of the power spectra, the data 
below the cut-off frequency play a minor role in determining a. Thus,a lower limit of 
/ is not included in the fitting algorithm. A comparison to the results obtained by 
increasing / d (upper triangle) shows a good agreement, except for the region where the 
exponent starts to deviate from its initial low value. In this region, the slightly smaller 
a for the case of increasing / d indicates hysteresis, which is in agreement with a former 
investigation l id . 

To check possible influence of the finite system size, I vary the particle number 
from 180 to 406, which corresponds to a range of area fraction from 0.277 to 0.624. As 
shown in Fig. [4](b), similar behavior of the power law exponent is found for all the area 
fractions explored. This demonstrates the robustness of using 1// noise as an order 
parameter for the melting transition. More specifically, the trend of a rapid deviation 
from the initially low a ~ —1.5 into a saturated region persists. Hence, I fit the data 
with a constant value followed by a parabolic growth, and determine the threshold f t h as 
the intersection point that minimizes the standard error. Quantitatively, the threshold 
frequency stays at ~ 1.3 Hz, with larger error bars for small cluster sizes. Note that the 
smallest cluster with N = 180 corresponds to a radius of roughly 7 particle diameters. 
For such a small “crystal”, the energy obtained from the container wall presumably 
leads to relatively strong fluctuations of the internal structures, and consequently large 
data scattering for the initial crystalline and amorphous states. For even smaller cluster 
size, the reduced probability for all particles to stay in one cluster hinders the statistical 
measurements performed here. 

In the following part of this section, I try to provide a clue to the emerging 1 // 
noise and its connections to the dynamics of the wet granular clusters in the vicinity of 
the melting transition. As a starting point, we need to identify the essential ingredients 
leading to the 1// noise. As already discussed in the above analysis, each assembly of 
the wet granular particles can be considered as one metastable state. This is reminiscent 
to the packing of equal sized spherical particles in three dimensions: The face-centered 
cube packing (fee) with a density of 0.74 is rare (50]. From time to time, we obtain 
a distribution of metastable states with a packing density range from random close to 
loose packings [5H [52]. 

The impacts with the container wall lead to frequent switches of the cluster from 
one metastable state to another one emerge, which can be considered as individual 
activation-relaxation processes. For example, the collision of a cluster in the amorphous 
state with the container wall leads to a certain energy injection E in j oc (1 — e 2 ), where 
e is the effective COR. On one hand, leads to an activation of defects inside the 
cluster, a corresponding reduction of P&, and an increase of the total potential energy 
through breaking of liquid bridges. On the other hand, it effectively “heats” the cluster 
up, i.e. enhances the kinetic energy of granular particles therein. The latter effect 
leads to a healing of defects generated, because the structure with four fold rotational 
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symmetry along the dislocation defects is locally unstable (see Fig. [2] and corresponding 
discussions). Moreover, the additional kinetic energy obtained from the healing process 
of a large defect may in turn generate small defects. A cascade of such energy transfer 
processes toward smaller and smaller length scales eventually leads to a relaxation of 
the injected energy through the dissipative interactions between individual particles. 
Following the BSM model [15, [16j [17, [18], we may assume that the relaxation process 
after an impact has the following autocorrelation function 

C(t) = A ■ e~ 7 o , (2) 


where A and r 0 denote the activation energy and the relaxation time, respectively. Note 
that the two parameters may differ dramatically from one impact to another, because the 
relaxed state after one impact is not necessarily the one with the lowest potential energy. 
Instead, a distribution of various metastable states arises with each state corresponding 
to a certain configuration of the wet granular particles. 

As C(t) is an even function, one can then calculate the corresponding power 
spectrum for this process with 

roo 

Pi(f) = 2 / Ae T o cos(27r/r)dr, (3) 

Jo 

which yields the Lorentzian function 


P(f) 


2Atq 

1 + (2tt/t 0 ) 2 ' 


(4) 


Assuming individual impacts are independent with each other, the power spectrum P(f) 
for the structure fluctuations measured can be obtained as a superposition of Lorentzian 
processes 

(5) 


where the probability Pr(A, To) of having activation energy A and relaxation time To 
crucially determines the exponent of P(/). For example, the distribution Pr(A, To) oc 
1 /to leads to P(f) oc 1/f for the whole frequency range. In reality, a limit on the 
relaxation time, say To G [ti,T 2 ], will lead to a limit on the power law decay: The 
spectrum flattens below 1 /t 2 and steepens above 1/ti. Therefore, one can determine 
the intrinsic time scale of the relaxation process from the power spectra. 

Note that the relaxation process through a cascade of defects nucleation and 
healing events relies on a certain rigidity of the cluster. As melting starts at relatively 
large /<j, the spatial as well the associated time correlation vanishes, and the energy 
injection and dissipation occur locally in the molten layers close to the impact point. 
As a consequence, a white noise type fluctuations with vanishing time correlation 
arise. As the internal structure fluctuations can be considered as a superposition of 
individual activation-relaxation processes, the dependency of the power law exponent 
on the stationary states of the nonequilibrium system represents the change of the 
corresponding correlation time. 
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So far, we have a qualitative understanding on why the power law exponent can be 
used to determine the transitions between various stationary states in the nonequilibrium 
model system. However, one has to face the following obstacles toward a quantitative 
understanding: There exists possible coupling between individual activation-relaxation 
events. A hard impact with the container wall will change the structure of a cluster 
substantially, which in turn will influence the COR for the next impact. Hence the 
assumption of independent subsequent events may not be justified. Moreover, such 
an influence will also lead to a certain distribution of the activation factor A, i.e. the 
energy injection, and a coupling between A and the relaxation time To- Therefore it is 
necessary to find the relation between A and To for more quantitative comparisons to 
the experiments. 

5. Conclusions 

In conclusion, I investigate the self-organization of a two-dimensional wet granular 
cluster under horizontal swirling motion, and use the percentage of particles in a local 
hexagonal packing to characterize the internal structure fluctuations with time. In the 
frequency domain, the power spectra of all fluctuations exhibit predominately a power 
law behavior with decay exponent ranging from a ~ —1.5 to ~ —0.3. More specifically, 
the exponent is closely related to the nonequilibrium stationary states of the system: For 
the initial crystalline and amorphous states at low driving frequencies, it stays low at 
~ —1.5, suggesting a diffusion noise. For the molten state at large driving frequencies, it 
saturates at ~ —0.3, suggesting a trend toward white noise. In the intermediate states 
associated with the melting transition, 1// type of noise with a ~ — 1 arises. A variation 
of the cluster size yields qualitatively the same power law behavior and quantitatively 
the same dependency of the decay exponent with /a. Such a robust behavior suggests an 
alternative way of detecting phase transitions in systems driven far from thermodynamic 
equilibrium. Finally, I show that the 1/f noise can be qualitatively understood from 
the impact mechanisms between the wet granular cluster and the container wall. 

This investigation demonstrates the importance of noise in characterizing systems 
driven out of thermodynamic equilibrium. Concerning the robustness of such a 
characterization, future questions to ask are: Will any type of fluctuation related to 
the energy flux in and out of the nonequilibrium system, e.g. the sound energy arising 
from the inelastic collisions, exhibits a similar behavior? Will the way of energy injection 
or dimensionality of the system influence the behavior? Concerning a more quantitative 
understanding of the emerging 1/f noise, it is also important to characterize the 
coupling between the activation energy and the relaxation time experimentally through 
a measurement of the effective coefficient of restitution for wet granular clusters. 
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